/*******************************************************************************

                                HKLREAD.C

  Title:        Hkl read package for Speden.
  Author:       Hanna Szoke
  Date:         5/18/2000 separated from hklutil.c
  Function:     This package contains functions that read structure factors:

			readfobs()			 
			readfcalc()
			readFcList()

		and low-level functions called by them.

*******************************************************************************/

#include	"util.h"
#include	"cellparams.h"
#include	"symmetry.h"
#include	"dims.h"
#include	"sp.h"

#define	FEPSI	1.e-4		/* for identifying close fractional (h,k,l) */
#define	DEPS	1.e-2 		/* criterion for judging relative errors in 
				structure factors that should agree */
#define	MAX_MISMATCH	20 	/* # of mismatches reported in non-verbose 
				mode */

real	amp0 = -1 ;			/* in dims.h 	*/
real	phase0 = -1 ;
real    F000 = 0 ;
real    S000 = 0 ;


struct Ntypes {		/* for reading exp. and calculated structure factors */
     int used;
     int total ;
} ;

static	void	makeSigmas() ;
static	int     set_off_grid_fc(real *, real *, float *, float *, float *) ;
static	int    	set_fo(int, int *, real, real, float, float, float) ;


	/*********************************************************
	Oversee reading of experimental (.fo) ungridded diffraction 
	pattern having format
		IND  %f  %f  %f  FOBS  %g  SIG  %g
	*********************************************************/
					
int readfobs(char *name)
{
     FILE	*fp ;
     int	first3chars(char *, char *) ;
     real	amplitude ;
     real	sigvalue ;
     struct	Ntypes	N ;
     int	Nlen ;
     int        n, m ;
     float	h, k, l ;
     char	generateSigmas = FALSE ;

     if (strcmp(name, "empty") == 0) 
	  Fatal("You cannot use 'empty' for Fobs file in SpEden") ;

     check_exist(name) ;
     Nlen = filelength(name) ;
     printf("Filelength = %d\n", Nlen) ;
 
     for (n = 0; n < Nlen; n++) {
          *(Folist+n) = 0. ;
          *(Siglist+n) = 0. ;
     }

     m = N.used = N.total = 0 ;

     fp = fopen(name, "r") ;

	/***********************************************************
	Read experimental data, expand to P1.
        Generate (h,k,l) triplets such that the Bragg condition
        is satisfied and at most only one of the two reciprocal
        lattice points -- (h,k,l) and (-h,-k,-l) -- appears.  
	Arbitrarily, we select h >= 0.  We use Friedel's law.
	***********************************************************/

     while (fgets(nextline, MAXSTRING, fp) != NULL)  {
	  if (((int)strlen(nextline) > 1) &&
	      (first3chars(nextline, "IND") == 0)) {

	     N.total++ ;
             sscanf(nextline, "%*s %f %f %f %*s %g %*s %g", 
		    &h, &k, &l, &t1, &t2) ;
	     amplitude = t1 ;
	     sigvalue = t2 ;  
             if (amplitude < 0) {
                   Warning(
	   "Your input fobs file contains entries with amplitude < 0;") ;
                   Fatal("Quitting.") ;
	      }
              if ((strstr(nextline, "SIG") != NULL)||
                  (strstr(nextline, "sig") != NULL))
		   generateSigmas = TRUE ;
              else
                  sigvalue = 0 ;

             if (m == 0) {	/* this should be where F000 & S000 are read */
                if ((h != 0) || (k != 0) || (l != 0)) {
                   Warning(
                 "Your first Fobs entry should be index (0 0 0)") ;
                   Fatal("Quitting.") ;
	        }
                else {
                   F000 = t1;
                   S000 = t2;
                }
             }
	     if (amplitude >= 0) {

                 N.used += set_fo(N.total, &m, amplitude, sigvalue, h, k, l) ;

             }		/* end check for potentially valid input */
         }		/* end check for something on input line */
     }			/* end check for EOF	*/

     fclose(fp) ;

     sprintf(message, "%d out of %d Fobs entries used;", N.used, N.total) ;
     printTwice(message) ;

     Nfob = N.used ;
     if (generateSigmas) {
         printTwice("Generating  sigmas") ;
         makeSigmas() ;
     }

     sprintf(message, "Fobs(0,0,0) = %g, sigma = %g", F000, S000) ;
     printTwice(message) ;
     
     return(N.used) ;
}
int readfobs_chunk(char *name, int Nfrom, int Nto)
{
     static	FILE	*fp ;
     int	first3chars(char *, char *) ;
     real	amplitude ;
     real	sigvalue ;
     struct	Ntypes	N ;
     int	Nlen ;
     int        n, m ;
     float	h, k, l ;
     char	generateSigmas = FALSE ;

     if (Nfrom == Nto) {
          fclose(fp) ;
          return(0) ;
     }

     if (strcmp(name, "empty") == 0) 
	  Fatal("You cannot use 'empty' for Fobs file in SpEden") ;

     check_exist(name) ;
     Nlen = Nto - Nfrom ;
     printf("Reading %d lines\n", Nlen) ;
 
     for (n = 0; n < Nlen; n++) {
          *(Folist+n) = 0. ;
          *(Siglist+n) = 0. ;
     }

     m = N.used = N.total = 0 ;

     if (Nfrom == 0) fp = fopen(name, "r") ;

	/***********************************************************
	Read experimental data, expand to P1.
        Generate (h,k,l) triplets such that the Bragg condition
        is satisfied and at most only one of the two reciprocal
        lattice points -- (h,k,l) and (-h,-k,-l) -- appears.  
	Arbitrarily, we select h >= 0.  We use Friedel's law.
	***********************************************************/

     while (fgets(nextline, MAXSTRING, fp) != NULL)  {
	  if (((int)strlen(nextline) > 1) &&
	      (first3chars(nextline, "IND") == 0)) {

	     N.total++ ;
             sscanf(nextline, "%*s %f %f %f %*s %g %*s %g", 
		    &h, &k, &l, &t1, &t2) ;
	     amplitude = t1 ;
	     sigvalue = t2 ;  
             if (amplitude < 0) {
                   Warning(
	   "Your input fobs file contains entries with amplitude < 0;") ;
                   Fatal("Quitting.") ;
	      }
              if ((strstr(nextline, "SIG") != NULL)||
                  (strstr(nextline, "sig") != NULL))
		   generateSigmas = TRUE ;
              else
                  sigvalue = 0 ;

             if ((m == 0) &&
                 (Nfrom == 0)){	/* this should be where F000 & S000 are read */
                if ((h != 0) || (k != 0) || (l != 0)) {
                   Warning(
                 "Your first Fobs entry should be index (0 0 0)") ;
                   Fatal("Quitting.") ;
	        }
                else {
                   F000 = t1;
                   S000 = t2;

                   sprintf(message, 
                           "Fobs(0,0,0) = %g, sigma = %g", F000, S000) ;
                   printTwice(message) ;
                }
             }
	     if (amplitude >= 0) {

                 N.used += set_fo(N.total, &m, amplitude, sigvalue, h, k, l) ;

             }		/* end check for potentially valid input */
         }		/* end check for something on input line */
         if (N.total == Nlen) 
             break ;
     }			/* end check for EOF	*/

     sprintf(message, "%d out of %d Fobs entries used;", N.used, N.total) ;
     printTwice(message) ;

     if (generateSigmas) {
         printTwice("Generating  sigmas") ;
         makeSigmas() ;
     }
     
     return(N.used) ;
}

	/*********************************************************
	Allocate space for fc list file.
	*********************************************************/

void alloc_fc(int N)
{
     strcpy(message, "alloc_fc") ;

     Fc_relist = (real *) e_malloc(N*sizeof(real), message) ;
     Fc_imlist = (real *) e_malloc(N*sizeof(real), message) ;
     Fc_Hlist =  (float *) e_malloc(N*sizeof(float), message) ;
     Fc_Klist =  (float *) e_malloc(N*sizeof(float), message) ;
     Fc_Llist =  (float *) e_malloc(N*sizeof(float), message) ;

     return ;
}

	/*********************************************************
	Oversee reading and generation of .fc ungridded 
	diffraction patterns.  
	*********************************************************/
					
int readFcList(char *name)
{
     FILE	*fp ;
     int	first3chars(char *, char *) ;
     real	amplitude ;
     real	phase ;
     struct	Ntypes	N ;
     int        n, m ;
     float	h, k, l ;

     if (strcmp(name, "empty") == 0) 
	  Fatal("You cannot use 'empty' for list Fc file in SpEden") ;

     check_exist(name) ;

     fp = fopen(name, "r") ;

     /* Initialize arrays */

     for (n = 0; n < Nfc; n++) {
          *(Fc_relist+n) = 0. ;
          *(Fc_imlist+n) = 0. ;
     }

     m = N.used = N.total = 0 ;

	/***********************************************************
	Read experimental data, expand to P1.
        Generate (h,k,l) triplets such that the Bragg condition
        is satisfied and at most only one of the two reciprocal
        lattice points -- (h,k,l) and (-h,-k,-l) -- appears.  
	Arbitrarily, we select h >= 0.  We use Friedel's law.
	***********************************************************/

     while (fgets(nextline, MAXSTRING, fp) != NULL)  {
	  if (((int)strlen(nextline) > 1) &&
	      (first3chars(nextline, "IND") == 0)) {

	     N.total++ ;
             sscanf(nextline, "%*s %f %f %f %*s %g %g", 
		    &h, &k, &l, &t1, &t2) ;
	     amplitude = t1 ;
	     phase = t2 ;  
             if (amplitude < 0) {
                   Warning(
	   "Your input fcalc file contains entries with amplitude < 0;") ;
                   Fatal("Quitting.") ;
	      }

                 N.used += set_off_grid_fc(&amplitude, &phase, &h, &k, &l) ;

         }		/* end check for something on input line */
     }			/* end check for EOF	*/

     fclose(fp) ;

     sprintf(message, "%d out of %d Fcalc entries read in;", N.used, Nfc) ;
     printTwice(message) ;

     sprintf(message, "Fcalc(0,0,0) = (%g, %g)", amp0, phase0) ;
     printTwice(message) ;
     
     Nfc = N.used ; 

     return(Nfc) ;
}

	/***********************************************************
	Generate sigmas using S(H) = sqrt(Fav / C + F(H)) / 10.  
	Apr 1, 2003: using C = 10. in place of 1.
	Sept. 20, 2005: if there's an S000, don't overwrite it!
	***********************************************************/

void	makeSigmas()
{
     int	n, n0 ;
     real	fav = 0 ;

     for (n = 0; n < Nfob; n++)		
          fav += *(Folist + n) ;

     fav /= Nfob ;

     n0 =  (S000 > 0.0) ? 1 : 0 ;
         
     for (n = n0; n < Nfob; n++)
          *(Siglist+n) = sqrt(fav/10. + *(Folist+n)) / 10.;

     return ;
}

	/***********************************************************
	Identify a position in the h >= 0 half ellipsoid.  
	***********************************************************/

int    set_fo(int N, int *m, real amplitude, real sigvalue, 
	       float h, float k, float l)
{
     real     fo_nusq(int) ;

         if (h < 0) {
	     h = -h ;
	     k = -k ;
	     l = -l ;
         }

         if (h > Nh-1)
             return (0) ;
         if ((k < -Nk/2) || (k >= Nk/2))
             return (0) ;
         if ((l < -Nl/2) || (l >= Nl/2))
             return (0) ;

	 *(Hlist + *m) = h ;
	 *(Klist + *m) = k ;
	 *(Llist + *m) = l ;

         if (fo_nusq(*m) < limit_rat) {

	      *(Folist + *m) = amplitude ;
	      *(Siglist + *m) = sigvalue ;
	      (*m)++ ;
	      return(1) ;
         }
	 else
	      return(0) ;
}

	/***********************************************************
	Identify a position in the h >= 0 half ellipsoid.  
	***********************************************************/

int    set_off_grid_fc(real *amp, real *phase, float *h, float *k, float *l)
{
     static	int	m = 0 ;
     static	int	mtot = 0 ;

    mtot++ ;

/****
     if (*h < 0) {
          *h = -*h ;
          *k = -*k ;
          *l = -*l ;
          sign = -1 ;
     }
     if ((*h < 0) || (*h >= Nh/2))
             return (0) ;
****/

/* I cannot understand this chunk of code so i'll comment it out and use the old version*/

/*     if ((*h < -Nh/2) || (*h >= Nh/2))
             return (0) ;
     if ((*k < -Nk/2) || (*k >= Nk/2))
             return (0) ;
     if ((*l < -Nl/2) || (*l >= Nl/2))
             return (0) ;
*/

     if (*h > Nh-1)
             return (0) ;
     if ((*k < -Nk/2) || (*k >= Nk/2))
             return (0) ;
     if ((*l < -Nl/2) || (*l >= Nl/2))
             return (0) ;


     if ((*h==0) && (*k==0) && (*l==0)) {
         amp0 = *amp ;
         phase0 = *phase ;
         if (m != 0) { 
            sprintf (message, 
		 "F(000) is the %d-th entry in your fcalc input.", mtot) ;
            Warning(message) ;
         }
     }

      *(Fc_Hlist + m) = *h ;
      *(Fc_Klist + m) = *k ;
      *(Fc_Llist + m) = *l ;

      *(Fc_relist + m) = *amp * cos(DTOR* *phase) ;
      *(Fc_imlist + m) = *amp * sin(DTOR* *phase) ;
      m++ ;
     
      return (1) ;
}
